Non-equilibrium dynamics of the Bose-Hubbard model: A projection operator 

approach 
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We study the phase diagram and non-equilibrium dynamics, both subsequent to a sudden quench 
of the hopping amplitude J and during a ramp J(t) = Jt/r with ramp time t, of the Bose-Hubbard 
model at zero temperature using a projection operator formalism which allows us to incorporate the 
effects of quantum fluctuations beyond mean-field approximations in the strong coupling regime. 
Our formalism yields a phase diagram which provides a near exact match with quantum Monte 
Carlo results in three dimensions. We also compute the residual energy Q, the superfluid order 
parameter A(t), the equal-time order parameter correlation function C(t), and the wavefunction 
overlap F which yields the defect formation probability P during non-equilibrium dynamics of the 
model. We find that Q, F, and P do not exhibit the expected universal scaling. We explain this 
absence of universality and show that our results compare well with recent experiments. 
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Ultracold atoms in optical lattices provide us with an 
unique setup to study non-equilibrium quantum dynam- 
ics of closed quantum systems [l], Q. The theoretical 
study of such quantum dynamics has seen great progress 
in recent years [3[. Most of these theoretical works have 
either restricted themselves to the physics of integrable 
and/or one-dimensional (ID) models or concentrated on 
generic scaling behavior of physical observables for sud- 
den or slow dynamics through a quantum critical point 
(QCP) However, quantum dynamics of specific ex- 
perimentally realizable non-integrable models in higher 
spatial dimension d and strong coupling regime have not 
been studied extensively. The Bose-Hubbard model with 
on-site interaction strength U and nearest neighbor hop- 
ping amplitude J, which provides an accurate descrip- 
tion for ultracold bosons in an optical lattice, constitutes 
an example of such models [4]. Most of the studies on 
dynamics of this model have concentrated on numerics 
for d < 2 0], weak coupling regime Q, and mean- 
field description of quench dynamics in the strong cou- 
pling regime [8, 9]. Recent experiments [2] on higher di- 
mensional Bose-Hubbard models in the strong-coupling 
regime (U ^> J) clearly necessitate computation of dy- 
namical evolution of several quantities beyond the mean- 
field theory and for arbitrary ramp time r. To the best 
of our knowledge, such a study has not been carried out. 

In this work we present a theoretical formalism beyond 
mean-field theory, which enables us to investigate in a 
semi-analytic way, at equal footing, both the equilibrium 
phase diagram and the non-equilibrium dynamics of the 
Bose-Hubbard model in the strong coupling regime and 
at zero temperature. The central results of our work are 
the following. First, we compute the equilibrium phase 
diagram and demonstrate that it provides a near per- 
fect match to the corresponding quantum Monte Carlo 
(QMC) results [To[ in three dimensions (3D). Second, we 
apply our formalism to non-equilibrium dynamics of the 



model for a finite ramp J{t) = Jt/r from t = U to t = tf. 
We compute the residual energy Q and the wavefunction 
overlap F [i.e. the overlap between the system wavefunc- 
tion after the ramp and the corresponding ground state 
wavefunction with J = J(t/)] which also yields the de- 
fect formation probability P = 1 — F [3] as functions of 
r. We show that for slow ramps P reaches a plateau, 
showing absence of expected scaling behavior [3[. We 
qualitatively explain such an absence of universal scal- 
ing and relate it to the recent experimental observations 
of Ref. |2[. Finally, we show that our formalism allows 
us to address the time evolution of the bosons after a 
sudden quench from the Mott ( J = Ji) to the superfluid 
(J = Jf) phase through the tip of the Mott lobe. We 
compute the order parameter A(t) and the equal-time 
order parameter correlation function C(t) during such 
an evolution. We also compute F and Q for a sudden 
quench from the critical point ( Ji = J c ) to the superfluid 
phase and show that they agree to the finite ramp results 
in the limit of small r and do not exhibit universal scal- 



ing behavior 11]. We note that dynamical properties of 
the Bose-Hubbard model for d > 2 in the strongly cou- 
pled regime have not been addressed beyond mean-field 
theory so far; our semi-analytical results therefore con- 
stitute significant extension of our understanding of the 
dynamics of this model in the strong-coupling regime. 
The Hamiltonian of the Bose-Hubbard model is 

n = ^-J6^+^[-/in r +|ri r (n r -l)], (1) 

(rr'> r 

where b (n) is the boson annihilation (number) opera- 
tor living on the sites of a <i-dimensional hypercubic lat- 
tice, and the chemical potential \i fixes the total number 
of particles. The corresponding many-body Schrodinger 
equation ihdt\ip) = H\^) is difficult to handle even nu- 
merically due to the infinite dimensionality of the Hilbert 
space. A typical practice is to use the Gutzwiller ansatz 
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FIG. 1: (Color online) (a) Schematic representation of the 
Mott state with ft — 1. (b) Typical hopping process mediated 
via T°. (c) Hopping process mediated via P^TgP^ . Notice 
that the states in (c) become part of the low-energy manifold 
near the critical point, while those in the right side of (b) do 
not and are always at an energy U above the Mott state. 
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FIG. 2: (Color online) Phase diagram of the Bose-Hubbard 
model in 2D (a) and 3D (b). The blue dots and blue solid lines 
(black dashed line) indicate the phase diagram obtained by 
the projection operator (mean- field) method. The red squares 
indicate QMC data. 



= exp(iS)\iJj). We use a Gutzwiller ansatz 

i^) = IIE/» r) i")' 



(3) 



= Yl r ^2n c ^\ n ) an d solve for Cn J keeping a finite 
number of states n around the Mott occupation number 
n = n. This yields the standard mean-field results with 
Cn^ = c n for homogeneous phases of the model [l2| . 

To build in fluctuations over such a mean-field the- 
ory, we use a projection operator technique [l3j |. The 
key idea behind this approach is to introduce a projec- 
tion operator Pi = |n)(n| r x |n)(n| r /, which lives on 
the link £ between the two neighboring sites r and r / '. 
The effect of Pi is to project any state of the system 
to the manifold of states for which n r ,n r / = n. Us- 
ing Pi : one can rewrite the hopping term of H: T' — 
£< rr <) -Jbtbr< =T,e T t = EeKPeTt + TtPd + PfTePf], 
where P^ = (1 — Pi). Note that, as schematically ex- 
plained in Fig.[TJ in the strong-coupling regime, the term 
Tg[J] = {PgTi+TgPi) represents hopping processes which 
take the system out of the low-energy manifold [13|- To 
obtain an effective low energy Hamiltonian, we there- 
fore devise a canonical transformation via an operator 
S = S[J] = ^ £ i[Pi,Ti]/U, which eliminates T°[J] to 
first order in zqJ/U, where zq = 2d is the coordination 
number of the lattice. This leads to the effective Hamil- 
tonian H* = exp(iS)Hexp(-iS) up to 0(z$J 2 /U) 

t t 

J] H \PiTiTgt - TtPtTv 



PiTfPf - TpPpTp 



V) 



(2) 



where Ho denotes the on-site terms in Eq. HJ Us- 
ing H* one can now compute the ground state en- 
ergy E = (tl>\H\tl>) = (ip'\H*W) + 0(^J 3 //7 2 ), where 



so that = only in the Mott limit (S, J = 0) and 
the energy becomes a functional of the coefficients 



E[{fn}]J] = (4>'\H*\4>'). 



(4) 



The contributions to E[{f n }; J] from the first two terms 
in H* (Eq. [2]) represent the mean- field energy functional, 
while the rest of the terms yield contributions due to 
quantum fluctuations. Thus the method constitutes sys- 
tematic inclusion of effects of quantum fluctuation over 
mean-field theory. The phase diagram obtained by min- 
imizing E[{f n };J] with respect to {f n } for 2D(3D) and 
n = 1 is shown in Fig. [2ja)(Fig. [2(b)). We find that the 
match with QMC data [10| is nearly perfect in 3D (with 
an error of ~ 0.05% at the tip of the Mott lobe) where 
mean- field theory provides an accurate starting point. In 
contrast, for 2D, we find J c /U = 0.055 compared to the 
QMC value 0.061 (red line in Fig. 12(a)). Here the match 
with QMC is not as accurate; however it compares fa- 
vorably to other analytical methods [15| . For the rest of 
this work, we shall restrict ourselves to d = 3 and n = 1. 

Next, we apply our formalism to address the dynamics 
of the model during a ramp with finite rate r _1 . We con- 
sider a ramp process under which J evolves from Ji at 
ti = to Jf at tf = r:J(t) = J{ + (Jf — Ji)t/r. To solve 
the Schrodinger equation, we make a time-dependent 
canonical transformation via a time-dependent S[J(t)] 
to eliminate T® up to first order from H at each instant. 
This yields the Schrodinger equation 



(ihdt + dS/&t)\rl>') = H*[J(t)]W) 



(5) 



The additional term dS/dt takes into account the possi- 
bility of creation of excitations during the time evolution 
with a finite ramp rate r _1 . The above equation yields 
an accurate description of the ramp with H* [J(t)] given 
by Eq.[2]for J(t)/U <C 1. Note that this does not impose 
a constraint on r; it only restricts Jf/U and Ji/U, to 
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FIG. 3: (Color online) (a) Plot of P as a function rll (in units 
of h = 1) for Ji/U = 0.05 (SF phase) and Jf/U = 0.005 (Mott 
phase) showing the plateau-like behavior at large r, and the 
corresponding saturation of Q (b). The inset in (b) shows Q 
and 1 — F as a function of SJ/J C for rll = 1. 



be small. Thus the method can treat both "slow" and 
"fast" ramps at equal footing. 

Substituting Eq. [3] into Eq. [5] and allowing for time- 
dependent fn^ , we find that the evolution of the system 
is given by the set of coupled equations 



ihdtfP = SE[{f n (t)}; J(t)]/J/*W + (6) 

+ V^+T/n+l ^n^/,n-l-^,n-l^n 

where ^ r = (^|6 r |^) = £ n *>m = £ n V^TT/n* (r) /Si, 
and S nn ' is the Kronecker delta. 

Using Eq. [6j we solve for ffi = / n for transla- 
tionally invariant systems numerically keeping all states 
< n < 5 with n = 1. Using these, we compute the de- 
fect formation probability P = 1 — F = 1 — | (^g I ^ (fy ) ) 1 2 ? 
where (IV^/))) denotes the final ground state (state 
after the ramp), for a ramp from Ji/U = 0.05 (superfluid 
phase) to Jf/U = 0.005 (Mott phase) as a function of 
r. We find that P exhibits a plateau like behavior at 
large r and do not display universal scaling as expected 
from generic theories of slow dynamics of quantum sys- 
tems near critical point |3j- This seems to be in qualita- 
tive agreement with the recent experiments presented in 
Ref. [2], where ramp dynamics of ultracold bosons from 
superfluid to the Mott region has been experimentally 
studied. Indeed, it was found, via direct measurement 
of n per site, that P displays a plateau like behavior 
similar to Fig. [3fa) [the inset displays the saturation for 
longer r). In Fig. [3^b), we show the analogous satura- 
tion and lack of universal scaling of the residual energy 
Q = f\7~L[J f]\ip f) — Eo[Jf], where Eo[Jf] denotes the 
ground state energy at J = Jf as obtained by minimizing 
E[{fn};Jf] inEq.Sl 

Such a lack of universality in the dynamics can be qual- 
itatively understood from the absence of contribution of 
the critical (k = 0) modes. In the strong-coupling regime 
(J/U C 1), the system can access the k = modes after 
a time T, which can be roughly estimated as the time 
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FIG. 4: (Color online) Plot of A(t) (a) and C(t) (b) as a 
function of tU, for Jf = 1.02J C . (c) The time period T of the 
oscillations of A(t). (d) Same as in (a) for Jf = 3.51J C . We 
have set ft = 1 for all plots. 



taken by a boson to cover the linear system dimension 
L. For typical small J (U = 1) in the Mott phase and 
near the QCP, T ~ 0(LH/J) can be very large. Thus for 
t < T, the dynamics, governed by local physics, which 
is well captured by our method, do not display critical 
scaling behavior. We note that in realistic experimental 
setups in the deep Mott limit [2], T may easily exceed 
the system lifetime making observation of universal scal- 
ing behavior impossible in such setups. 

We now apply this method to address the dynam- 
ics of the model after a sudden quench 14], from Ji 
(Mott phase) to Jf (superfluid phase) through the tip 
of the Mott lobe, where the dynamical critical expo- 
nent z = 1. The time evolution of the order param- 
eter A r (t) = (^)M^)> = ^ (t)\b' Y W (t)) , where 
b' Y = exp(iS[Jf])b r exp(— iS[Jf]), can then be expressed 
in terms of fn^ as 



A r (t) = Mt) + J/U^™[\fi r) \ 2 -fn\\ 



+ (n+1) \A r) \ 2 ~A%\ 2 



</V,n-l 



$r,n-l Ifr'fi + $ rn $r,n-l 



(7) 



where $ rn = y/(n + l)(n + 2)/* (r) /n+2- Note that the 
first term in Eq. [7] represents the mean- field result. The 
role of quantum fluctuations in the evolution of A r (t) be- 
comes evident in computing the equal-time order parame- 
ter correlation function C r {t) = (^ f (t)\b f r b f r \^ f (t))-Al(t). 
To compute A r (t) and C r (£), we solve Eq. [6] numeri- 
cally for a translationally invariant system. The resul- 
tant plot of A r (t) = A(t) is shown in Fig. H(a)[(d)] 
for Ji = and J f /J c = 1.02(J//J C = 3.51). We find 
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FIG. 5: (Color online) Plot of 1 — F and Q as a function of 
the 5 J for S J/ J c <C 1. The lines correspond to fits yielding a 
power 1 - F(Q) ~ (£J) ri(r2) with n ~ 0.89 and r 2 ~ 1.9. 



that near the critical point, A(t) displays oscillations 
with a single characteristic frequency [8j, while away 
from the critical point (Jf/J c — 3.51), multiple fre- 
quencies are involved in its dynamics. The time period 
T (Fig. E{c)) of these oscillations near J c is found, as 
a consequence of critical slowing down, to have a di- 
vergence T ~ \J f - Ji\~ zv = 5J-°- 35±0 - 05 , leading to 
zv = 0.35 ± 0.05 for d = 3 [3]. Finally, in Fig. gfb) we 
plot C r (t) = C(£) as a function of t, for Jf = 1.02J C . 
We find that \C(t)/ A 2 (t)\ may be as large as 0.5 at the 
tip of the peaks of A(t), which shows strong quantum 
fluctuations near the QCP. 

Next, we compute the wavefunction overlap F = 



\^f\^ c )Y 



W 



f ' 



i5[J/] p -i5[J c ] 



|^)| 2 for a sudden 



quench starting at the QCP. Here | < 0/)(| < C )) denotes the 
ground state wavefunction for J = Jf(J c ). We also com- 
pute the residual energy Q = {^ c \H[Jf)\^c) — EgIJj], 
and in Fig. [5] we plot 1 — F and Q for the homogeneous 
case as a function of 5 J = \ Jf — J c |, for SJ/J C < 0.2. A 
numerical fit of these curves yields 1 — F ~ JJ 0,89 , and 
Q ~ JJ 1,90 , which disagrees with the universal scaling 
exponents (1 - F ~ 5J dv and Q - 5J^ d+z >) expected 



from sudden dynamics across a QCP with z = 1 [11]. 
Note that our results for the sudden quench match with 
those for the ramp dynamics at small r, shown in the 
inset of Fig. [3jb). In particular, the exponents obtained 
from the two cases are nearly identical, reflecting accu- 
rate reproduction of fast ramp dynamics in the sudden 
quench limit. 

Finally, we estimate the range of physical tempera- 
tures for which the zero temperature theory is accurate. 
For typical lattice depths in the Mott or critical regimes, 
U ~ 2 kHz ~ 200nK [1]. This yields, in 3D, a melting 
temperature T* ~ 0.2 £7 = 40nK for the Mott phase and a 
critical temperature T c ~ J c — 35nK for the SF phase 
at the Mott tip [161 ]. This requires the system tempera- 
ture to be a few nano-Kelvins (and <C T*, T c ) , which is 
well within the current experimental limit ~ InK [l6| . 



In conclusion, we have presented a projection operator 
formalism that describes in a semi-analytical way both 
the phase diagram, and non-equilibrium dynamics of the 
Bose-Hubbard model. It produces a phase diagram which 
is nearly identical to the QMC results in 3D, and allows 
computation of several quantities such as F, Q, A(t), 
P, and C(t) for non-equilibrium dynamics. Its predic- 
tion for P for a slow ramp matches qualitatively with 
recent experiments. The method, in principle, can be 
generalized to correlated systems which allow perturba- 
tive treatment of fluctuations and for studying ultracold 
bosons in a finite trap. We leave such considerations for 
future study. 
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